\(~\)
This notebook helps visualize and explore surface and groundwater observations stored in the National Water Information System (NWIS). It was written to perform automated pulls and interactive mapping and plotting of data as it was being collected and analyzed for a specific project, but it can be modified for any region where there is data in NWIS.
\(~\)
# watershed screening sites
screensites <- whatNWISsites(bBox = c(-70.448, 41.626, -70.392, 41.677))
nums <- "504|505|506|507"
screensites <- screensites %>%
filter(grepl(nums, station_nm))
# local demonstration sites
sites <- whatNWISsites(bBox = c(-70.405, 41.669, -70.395, 41.677))
# site info
baseurl1 <- "https://nwis.waterdata.usgs.gov/usa/nwis/qwdata/?site_no="
baseurl2 <- "&agency_cd=USGS"
siteinfo <- readNWISsite(sites$site_no) %>%
dplyr::select(site_no, station_nm, site_tp_cd, dec_lat_va, dec_long_va, dec_coord_datum_cd, alt_va, alt_datum_cd, well_depth_va)%>%
mutate(id = substr(station_nm, 9, 11)) %>%
mutate(weblink = paste(baseurl1,site_no,baseurl2,sep = " "))
# reclassify site types for mapping
siteinfo$site_tp_cd[grep("M01", siteinfo$station_nm)] <- "multilevel sampler (MLS)"
clust <- "505-0|517-0|524-0|525-0|526-0"
siteinfo$site_tp_cd[grepl(clust, siteinfo$station_nm)] <- "well cluster"
siteinfo$site_tp_cd[which(siteinfo$site_tp_cd=="GW")] <- "water table well"
siteinfo$site_tp_cd[which(siteinfo$site_tp_cd=="LK")] <- "surface water (pond)"
# add a treatment variable
siteinfo$treatment <- NA
controls <- c(505,508,514,515,520,521)
siteinfo$treatment[which(siteinfo$id %in% controls)] <- "control"
impacts <- c(517,522,523,524,525,526)
siteinfo$treatment[which(siteinfo$id %in% impacts)] <- "impact"
siteinfo$treatment[siteinfo$id == 509] <- "other"
# map all sites by type
pal <- colorFactor("Accent", domain = siteinfo$site_tp_cd)
leaflet(siteinfo, options = leafletOptions(zoomControl = FALSE)) %>%
addProviderTiles('Esri.WorldTopoMap')%>%
addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1,
fillOpacity = 1, fillColor = ~pal(site_tp_cd),
popup = ~paste(station_nm,"<br>","<a href='", weblink, "' target='_blank'>", "Link to data</a>")) %>%
addScaleBar("bottomright") %>%
addLegend(pal = pal, title = "Site type", opacity = 1, values = ~site_tp_cd, position = "topleft")
# map all sites by treatment
pal <- colorFactor("Accent", domain = siteinfo$treatment)
leaflet(siteinfo, options = leafletOptions(zoomControl = FALSE)) %>%
addProviderTiles('Esri.WorldTopoMap')%>%
addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1,
fillOpacity = 1, fillColor = ~pal(treatment),
popup = ~paste(station_nm,"<br>","<a href='", weblink, "' target='_blank'>", "Link to data</a>")) %>%
addScaleBar("bottomright") %>%
addLegend(pal = pal, title = "Treatment", opacity = 1, values = ~treatment, position = "topleft")
\(~\)
# parameters (nutrients samples are all filtered)
temp <- "00010" # temperature (degC)
sc <- "00095" # SpC (uS/cm)
ph <- "00400" # pH
do <- "00300" # DO (mg/L)
nh34 <- "00608" # ammonia and ammonium, as N (mg/L)
no2 <- "00613" # nitrite as N (mg/L)
no3 <- "00618" # nitrate as N (mg/L)
no23 <- "00631" # nitrate + nitrite as N (mg/L)
tn <- "62854" # total nitrogen (mg/L)
po4 <- "00671" # orthophosphate as P (mg/L)
delH <- "82082" # delta H2/H1 (per mil)
delO <- "82085" # delta O18/O16 (per mil)
pset <- c(temp, sc, ph, do, nh34, no2, no3, no23, tn, po4, delH, delO)
# RETRIEVAL
# pull water quality data
## option A (to be retired)
data <- readNWISqw(siteNumber = siteinfo$site_no, parameterCd = pset)
b <- data %>% dplyr::select(c(site_no, sample_dt, sample_tm, parm_cd, result_va))
# ## option B (suggested going forward)
# data <- readWQPqw(siteNumber = paste0("USGS-", siteinfo$site_no), parameterCd = pset)
# # select columns and rename
# b <- data %>% dplyr::select(c(MonitoringLocationIdentifier, ActivityStartDate, USGSPCode, CharacteristicName, ResultMeasure.MeasureUnitCode, ResultMeasureValue))
# names(b) <- c("site_no", "sample_dt", "parm_cd", "charname", "units", "result_va")
# b$site_no <- sub("USGS-", "", b$site_no)
# combine site info and qw data
a <- siteinfo
a$station_nm[grep("SHUBAEL", a$station_nm)] <- "SHUBAEL POND"
alldata <- right_join(a,b) %>%
rename(date = sample_dt) %>%
mutate(param = parm_cd, yearmon = format(date, "%Y-%m"))
# create param names from codes
alldata$param <- dplyr::recode(alldata$param, `00010` = "temp", `00095` = "SpC", `00400` = "pH", `00300` = "DO", `00608` = "NH34", `00613` = "NO2", `00618` = "NO3", `00631` = "NO23", `62854` = "TN", `00671` = "PO4", `82082` = "delH", `82085` = "delO")
# pull water levels data
gwl <- readNWISgwl(siteNumbers = siteinfo$site_no)
# join site info
dtw <- gwl %>%
filter(!is.na(lev_va))
dtw <- right_join(siteinfo %>% dplyr::select(station_nm, site_no, id, dec_lat_va, dec_long_va, alt_va, well_depth_va, weblink), dtw)%>%
rename(date = lev_dt)
# SUMMARIES
# number of sample depths at each location
nwells <- alldata %>%
group_by(id) %>%
summarise(depths = length(unique(well_depth_va)))
# number of depths and samples at each location, by parameter
allsum <- alldata %>%
group_by(id, param) %>%
summarise(count = n()) %>%
pivot_wider(id_cols = id, names_from = param, values_from = count)
allsum <- left_join(allsum, nwells) %>%
dplyr::select(id, depths, temp, SpC, pH, DO, PO4, NH34:NO3, TN) %>%
arrange(depths)
# estimate depth below water table (satdepth) for each sample
alldata <- left_join(alldata, dtw %>% select(c(date, site_no, lev_va))) # first, join depth to water
mls <- alldata %>% filter(site_tp_cd == "multilevel sampler (MLS)") %>% # for mls, assign depth of most shallow sampled port, -0.5ft, as estimate of depth to water
group_by(id) %>%
summarise(dtw = min(well_depth_va)-0.5)
for(i in 1:nrow(mls)){alldata$lev_va[alldata$id==mls$id[i]] = mls$dtw[i]}
alldata <- alldata %>%
mutate(satdepth = well_depth_va - lev_va)
# subset data for a specific sampling event
mapdate <- "2021-11-01"
zt <- alldata %>%
dplyr::filter(date > mapdate)
# number of depths and samples at each location, by parameter
ztsum <- zt %>%
group_by(id, param) %>%
summarise(count = n()) %>%
pivot_wider(id_cols = id, names_from = param, values_from = count)
ztsum <- left_join(ztsum, nwells) %>%
dplyr::select(id, depths, temp, SpC, pH, DO, PO4, NH34:NO3, TN) %>%
arrange(depths)
# means at each location, single time
zmean <- zt %>%
filter(!station_nm=="MA-A1W 517-0036") %>% # remove duplicate shallow well at A517
#filter(satdepth < 10) %>%
group_by(id, parm_cd) %>%
summarise(station_nm=first(id), dec_lat_va=median(dec_lat_va), dec_long_va=median(dec_long_va), weblink=first(weblink), well_depth_va=mean(well_depth_va, na.rm=TRUE), result_va=round(mean(result_va, na.rm=TRUE), 1))
# means at each location, time series
tsmean <- alldata %>%
filter(!station_nm=="MA-A1W 517-0036") %>% # remove duplicate shallow well at A517
#filter(satdepth < 10) %>%
group_by(id, yearmon, param) %>%
summarise(station_nm=first(id), dec_lat_va=median(dec_lat_va), dec_long_va=median(dec_long_va), weblink=first(weblink), well_depth_va=mean(well_depth_va, na.rm=TRUE), result_va=round(mean(result_va, na.rm=TRUE), 1), treatment=first(treatment))
# max over sample depths at each location
zmax <- zt %>%
group_by(id, parm_cd) %>%
summarise(station_nm=first(id), dec_lat_va=median(dec_lat_va), dec_long_va=median(dec_long_va), weblink=first(weblink), well_depth_va=mean(well_depth_va, na.rm=TRUE), result_va=max(result_va, na.rm=TRUE))
# calculate percentages of N in different forms
ftn <- pivot_wider(zt, id_cols = c(id, site_no:dec_long_va, well_depth_va), names_from = param, values_from = result_va)%>%
filter(!is.na(TN)) %>%
mutate(percent_organicN = round(100*(1 - ((NH34+NO23)/TN)), 1), percent_NO23 = round(100*NO23/TN, 1))
kable(ftn %>% arrange(id) %>% dplyr::select(station_nm, site_tp_cd, temp, SpC, pH, DO, PO4, NH34, NO3, NO23, TN, percent_organicN, percent_NO23), digits = 2)
| station_nm | site_tp_cd | temp | SpC | pH | DO | PO4 | NH34 | NO3 | NO23 | TN | percent_organicN | percent_NO23 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| MA-A1W 505-0036 | well cluster | 12.2 | 143 | 5.1 | 9.30 | NA | 0.02 | 3.18 | 3.18 | 3.23 | 0.9 | 98.5 |
| MA-A1W 508-M01-06WT | multilevel sampler (MLS) | 13.7 | 105 | 5.0 | 3.00 | NA | 0.02 | 3.28 | 3.28 | 3.38 | 2.4 | 97.0 |
| MA-A1W 509-0019 | water table well | 16.8 | 105 | 6.3 | 0.50 | 0 | 0.02 | 0.27 | 0.28 | 0.35 | 14.0 | 80.3 |
| MA-A1W 514-0028 | water table well | 11.7 | 138 | 4.7 | 5.50 | 0 | 0.02 | 6.14 | 6.14 | 6.30 | 2.2 | 97.5 |
| MA-A1W 515-0036 | water table well | 11.2 | 214 | 4.9 | 7.80 | NA | 0.02 | 11.90 | 11.90 | 12.40 | 3.9 | 96.0 |
| MA-A1W 517-0037 | well cluster | 11.3 | 134 | 4.8 | 5.50 | NA | 0.15 | 3.65 | 3.65 | 3.97 | 4.3 | 91.9 |
| MA-A1W 520-0048 | water table well | 10.7 | 134 | 5.3 | 7.80 | NA | 0.02 | 6.56 | 6.56 | 6.86 | 4.1 | 95.6 |
| MA-A1W 521-0041 | water table well | 10.7 | 133 | 5.2 | 9.60 | NA | 0.02 | 1.58 | 1.58 | 1.68 | 4.8 | 94.0 |
| MA-A1W 522-M01-06WT | multilevel sampler (MLS) | 16.4 | 65 | 5.5 | 0.23 | 0 | 0.02 | 0.18 | 0.18 | 0.23 | 12.6 | 78.7 |
| MA-A1W 523-M01-06WT | multilevel sampler (MLS) | 15.4 | 260 | 4.5 | 0.09 | 0 | 1.73 | 14.10 | 14.10 | 16.70 | 5.2 | 84.4 |
| MA-A1W 524-0039 | well cluster | 10.5 | 151 | 5.0 | 3.10 | NA | 0.02 | 3.05 | 3.05 | 3.24 | 5.2 | 94.1 |
| MA-A1W 525-0042 | well cluster | 11.7 | 183 | 4.8 | 9.10 | NA | 0.02 | 8.50 | 8.50 | 8.66 | 1.6 | 98.2 |
| MA-A1W 526-0034 | well cluster | 11.9 | 133 | 5.8 | 0.67 | NA | 0.02 | 0.42 | 0.43 | 0.53 | 15.8 | 80.4 |
\(~\)
\(~\)
\(~\)
# function to make the maps
mapparam <- function(data, palette, title){
leaflet(data, options = leafletOptions(zoomControl = FALSE)) %>%
addProviderTiles('Esri.WorldTopoMap')%>%
addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1,
fillOpacity = 1, fillColor = ~pal(result_va),
popup = ~paste(station_nm,"<br>","Value",result_va,"<br>","<a href='"
, weblink
, "' target='_blank'>"
, "Link to data</a>")) %>%
addScaleBar("bottomright") %>%
addLegend(pal = palette, title = title, opacity = 1, values = ~result_va, position = "topleft")
}
z <- zmean
ttl <- "Temp (°C)"
mapdata <- z %>% filter(parm_cd == temp)
pal <- colorNumeric("viridis", domain = mapdata$result_va)
m1 <- mapparam(mapdata, pal, ttl)
#m1
ttl <- "SpC (uS/cm)"
mapdata <- z %>% filter(parm_cd == sc)
pal <- colorBin("viridis", domain = mapdata$result_va)
m2 <- mapparam(mapdata, pal, ttl)
#m2
ttl = "pH"
mapdata <- z %>% filter(parm_cd == ph)
pal <- colorBin("viridis", domain = mapdata$result_va)
m3 <- mapparam(mapdata, pal, ttl)
#m3
ttl = "DO (mg/L)"
mapdata <- z %>% filter(parm_cd == do)
pal <- colorBin("viridis", domain = mapdata$result_va)
m4 <- mapparam(mapdata, pal, ttl)
#m4
ttl = "Nitrate, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == no3)
pal <- colorBin("viridis", domain = mapdata$result_va)
m5 <- mapparam(mapdata, pal, ttl)
#m5
ttl = "NH34, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == nh34)
pal <- colorBin("viridis", domain = mapdata$result_va)
m6 <- mapparam(mapdata, pal, ttl)
#m6
ttl = "Total N (mg/L)"
mapdata <- z %>% filter(parm_cd == tn)
pal <- colorBin("viridis", domain = mapdata$result_va)
m7 <- mapparam(mapdata, pal, ttl)
#m7
ttl = "Phosphate, as P <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == po4)
pal <- colorBin("viridis", domain = mapdata$result_va)
m8 <- mapparam(mapdata, pal, ttl)
#m8
# map depth to water
pal <- colorNumeric("viridis", domain = dtw$lev_va)
m9 <- leaflet(dtw %>% filter(date > mapdate), options = leafletOptions(zoomControl = FALSE)) %>%
addProviderTiles('Esri.WorldTopoMap')%>%
addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1,
fillOpacity = 1, fillColor = ~pal(lev_va),
popup = ~paste(station_nm,"<br>","Value",lev_va,"<br>","<a href='", weblink, "' target='_blank'>", "Link to data</a>")) %>%
addScaleBar("bottomright") %>%
addLegend(pal = pal, title = "Depth to water", opacity = 1, values = ~lev_va, position = "topleft")
sync(m1, m2, m3, m4, m5, m6, m7, m8, m9, ncol = 3, sync.cursor = FALSE)
\(~\)
\(~\)
\(~\)
z <- zmax
ttl <- "Temp (°C)"
mapdata <- z %>% filter(parm_cd == temp)
pal <- colorNumeric("viridis", domain = mapdata$result_va)
m1 <- mapparam(mapdata, pal, ttl)
#m1
ttl <- "SpC (uS/cm)"
mapdata <- z %>% filter(parm_cd == sc)
pal <- colorBin("viridis", domain = mapdata$result_va)
m2 <- mapparam(mapdata, pal, ttl)
#m2
ttl = "pH"
mapdata <- z %>% filter(parm_cd == ph)
pal <- colorBin("viridis", domain = mapdata$result_va)
m3 <- mapparam(mapdata, pal, ttl)
#m3
ttl = "DO (mg/L)"
mapdata <- z %>% filter(parm_cd == do)
pal <- colorBin("viridis", domain = mapdata$result_va)
m4 <- mapparam(mapdata, pal, ttl)
#m4
ttl = "Nitrate, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == no3)
pal <- colorBin("viridis", domain = mapdata$result_va)
m5 <- mapparam(mapdata, pal, ttl)
#m5
ttl = "NH34, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == nh34)
pal <- colorBin("viridis", domain = mapdata$result_va)
m6 <- mapparam(mapdata, pal, ttl)
#m6
ttl = "Total N (mg/L)"
mapdata <- z %>% filter(parm_cd == tn)
pal <- colorBin("viridis", domain = mapdata$result_va)
m7 <- mapparam(mapdata, pal, ttl)
#m7
ttl = "Phosphate, as P <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == po4)
pal <- colorBin("viridis", domain = mapdata$result_va)
m8 <- mapparam(mapdata, pal, ttl)
#m8
sync(m1, m2, m3, m4, m5, m6, m7, m8, m9, ncol = 3, sync.cursor = FALSE)
\(~\)
\(~\)
\(~\)
ttl = "Nitrate, as N, (mg/L) in groundwater wells and pond"
p <- ggplot(alldata %>% filter(site_tp_cd %in% c("well cluster", "water table well", "surface water (pond)")) %>% filter(parm_cd==no3), aes(x = date, y = result_va, group = station_nm)) +
geom_point(aes(color=id)) +
geom_line(aes(color=id)) +
theme_bw()+
ggtitle(ttl)
ggplotly(p)
ttl = "Nitrate, as N, (mg/L) in multilevel groundwater samplers (MLS)"
p <- ggplot(alldata %>% filter(site_tp_cd %in% c("multilevel sampler (MLS)")) %>% filter(parm_cd==no3), aes(x = date, y = result_va, group = station_nm)) +
geom_point(aes(color=id)) +
geom_line(aes(color=id)) +
theme_bw()+
ggtitle(ttl)
ggplotly(p)
ttl = "Nitrate, as N, (mg/L) in wells and MLS ports"
p <- ggplot(alldata %>% filter(param=="NO3" & !is.na(treatment)), aes(x = date, y = result_va, group = station_nm)) +
geom_point(aes(color = treatment)) +
geom_line(aes(color = treatment)) +
theme_bw()+
ggtitle(ttl)
ggplotly(p)
ttl = "Mean nitrate, as N, (mg/L) at site"
p <- ggplot(tsmean %>% filter(param=="NO3" & !is.na(treatment)), aes(x = yearmon, y = result_va, group = station_nm)) +
geom_point(aes(color = treatment)) +
geom_line(aes(color = treatment)) +
theme_bw()+
ggtitle(ttl)
ggplotly(p)
\(~\)
\(~\)
profdate <- "2021-06-01"
parameters <- c( "NO3") #"DO", "pH", "delO", "temp", "SpC"
profiledata <- alldata %>%
filter(id %in% c(505, 508, 517, 522, 523, 524, 525, 526)) %>%
filter(param %in% parameters) %>%
filter(date > profdate) %>%
filter(!station_nm=="MA-A1W 517-0036") %>% # remove duplicate shallow well at A517
#mutate(date = as.character(date)) %>%
#mutate(result_va = ifelse(param == "SpC", result_va/10, result_va)) %>%
#mutate(param = replace(param, param == "SpC", "SpC/10")) %>%
arrange(id)
xr <- c(0, 50)
yr <- c(80, 0)
# function to make vertical profiles
zprofile <- function(data, title, xrange, yrange){
p <- ggplot(data %>% group_by(yearmon), aes(x = result_va, y = well_depth_va)) +
geom_path(aes(color = yearmon, linetype = yearmon))+
scale_x_continuous()+ #limits = xrange
theme_bw()+
scale_y_reverse()+ #limits = yrange
ggtitle(title)
ggplotly(p)
}
# iterated over sites
p <- htmltools::tagList()
for (i in unique(profiledata$id)){
a <- filter(profiledata, id == i)
if(i==508){a <- filter(a, site_tp_cd == "multilevel sampler (MLS)")}
ttl <- paste0(a$site_tp_cd[1], " at ", i)
p[[i]] <- zprofile(a, ttl, xr, yr)
}
p